# ===============================================
# Paris ratchet paper
# Code for non-row-standardized trade and IO matrices
# ===============================================

#### Table of contents ####

#' 1. Load imputed dyadic data to remake trade and IOs matrices
#' 2. Run trade models
#' 3. Run IO models
#' 4. Interact trade and IOs
#' 5. Regional dummies


#### Initialize #### 

library(tidyverse)
library(knitr)
library(DescTools) # Winsorize()
library(broom)
library(modelsummary)
library(marginaleffects)
# library(readxl)
# library(lubridate)
# library(data.table)
# library(Amelia) # Multiple imputation
# library(lsa) # for cosine similarity
# library(WDI)
# library(magrittr)

options(scipen=999) # Scientific notation off
# options(scipen=0)  # Scientific notation on


## Set working directory

## Define core ISO3C codes to use in analysis
core_iso3c <- c("AFG", "AGO", "ALB", "AND", "ARE", "ARG", "ARM", "ATG", "AUS", "AUT", "AZE", "BDI", "BEL", "BEN", "BFA", "BGD", "BGR", "BHR", "BHS", "BIH", "BLR", "BLZ", "BOL", "BRA", "BRB", "BRN", "BTN", "BWA", "CAF", "CAN", "CHE", "CHL", "CHN", "CIV", "CMR", "COD", "COG", "COK", "COL", "COM", "CPV", "CRI", "CUB", "CYP", "CZE", "DEU", "DJI", "DMA", "DNK", "DOM", "DZA", "ECU", "EGY", "ERI", "ESP", "EST", "ETH", "FIN", "FJI", "FRA", "FSM", "GAB", "GBR", "GEO", "GHA", "GIN", "GMB", "GNB", "GNQ", "GRC", "GRD", "GTM", "GUY", "HND", "HRV", "HTI", "HUN", "IDN", "IND", "IRL", "IRN", "IRQ", "ISL", "ISR", "ITA", "JAM", "JOR", "JPN", "KAZ", "KEN", "KGZ", "KHM", "KIR", "KNA", "KOR", "KWT", "LAO", "LBN", "LBR", "LBY", "LCA", "LIE", "LKA", "LSO", "LTU", "LUX", "LVA", "MAR", "MDA", "MDG", "MDV", "MEX", "MHL", "MKD", "MLI", "MLT", "MMR", "MNE", "MNG", "MOZ", "MRT", "MUS", "MWI", "MYS", "NAM", "NER", "NGA", "NIC", "NLD", "NOR", "NPL", "NRU", "NZL", "OMN", "PAK", "PAN", "PER", "PHL", "PLW", "PNG", "POL", "PRK", "PRT", "PRY", "QAT", "ROU", "RUS", "RWA", "SAU", "SDN", "SEN", "SGP", "SLB", "SLE", "SLV", "SOM", "SRB", "SSD", "STP", "SUR", "SVK", "SVN", "SWE", "SWZ", "SYC", "SYR", "TCD", "TGO", "THA", "TJK", "TKM", "TLS", "TON", "TTO", "TUN", "TUR", "TUV", "TZA", "UGA", "UKR", "URY", "USA", "UZB", "VCT", "VEN", "VNM", "VUT", "WSM", "YEM", "ZAF", "ZMB", "ZWE")

## Define a "max" function that works with NA
my.max <- function(x) ifelse( !all(is.na(x)), max(x, na.rm=T), NA)



#### Make non--row-standardized spatial weight for trade ####

#### Load processed data from "00_make_analysis_data.R" to remake spatial weights

## Load panel of country covariates
country_level_average <- readRDS("output/country_level_average_20230920.rds")

## Load the dyadic imputed dataset --- trade variables
full_imputed_averaged <- readRDS("output/full_imputed_averaged_20230920.rds")

## Load main analysis dataset --- GDP 
data_for_models <- readRDS("output/data_for_models20230920.Rds")

## Partners targets
partners_targets_cs <- full_imputed_averaged %>% 
  select(iso3_o, iso3_d, year, income_d,
         # should all be _d's targets --- 6? + safe_
         # -- Percentage
         ndc1_percentage_average_d,
         safe_ndc1_percentage_average_d,
         # -- Laws
         climate_laws_d,
         climate_laws_new20142019_d) %>% 
  filter(iso3_o!=1,
         iso3_d!=1) %>% 
  filter(year %in% c(2015:2019)) %>% 
  group_by(iso3_d) %>%
  mutate(climate_laws_d = mean(climate_laws_d)) %>% 
  distinct() %>% 
  select(-iso3_o, -income_d) %>% 
  group_by(iso3_d) %>%
  summarise_at(vars(ndc1_percentage_average_d:climate_laws_new20142019_d),
               ~mean(.))


#' Then redo the process from "00_make_analysis_data.R" to make the spatial weights, 
#' but don't row-standardize
hold_dyadic <- full_imputed_averaged %>% 
  select(iso3_o, iso3_d, year, 
         tradeflow_baci, trade_eite) %>% 
  filter(iso3_o!=1,
         iso3_d!=1) %>% 
  # Summarize dyadic flows in the five year period
  filter(year>2015, year < 2020) %>% 
  group_by(iso3_o, iso3_d) %>%
  summarize(tradeflow_baci = mean(tradeflow_baci), # <- these are exports from _o to _d
            trade_eite = mean(trade_eite)) %>% 
  # For each dyad, get their total trade
  rowwise() %>% 
  mutate(dyad_id = paste(min(iso3_d, iso3_o),
                         max(iso3_d, iso3_o),
                         sep = "_")) %>% 
  group_by(dyad_id) %>% 
  mutate(dyadic_tradeflow = sum(tradeflow_baci), # Sum of exports/imports (undirected) in that dyad
         dyadic_trade_eite = sum(trade_eite)) %>% 
  group_by(iso3_o) %>% 
  mutate(exports_o_tradeflow = sum(tradeflow_baci), # Sum of _o's exports
         exports_o_eite = sum(trade_eite)) %>% 
  group_by(iso3_d) %>% 
  mutate(imports_d_tradeflow = sum(tradeflow_baci), # Sum of _d's imports; this is at the wrong level, re-shape below
         imports_d_eite = sum(trade_eite)) %>% 
  arrange(dyad_id) %>% 
  ungroup() %>% 
  as.data.frame()

## Second step of the transformation 
imports_iso3_o <- hold_dyadic %>% 
  group_by(iso3_d) %>% # Yes, should be _d; this is for re-shaping the above
  # Get country-level total imports in each category
  distinct(imports_d_tradeflow, 
           imports_d_eite) %>% # Yes, should be _d
  rename(iso3_o = iso3_d, # Then, rename and merge back in
         imports_o_tradeflow = imports_d_tradeflow,
         imports_o_eite = imports_d_eite) %>% 
  ungroup() %>% 
  as.data.frame()


## Extract GDP from country_level_average 
gdp_levels <- country_level_average %>% 
  filter(year == 2016) %>%
  # Imports and exports are already log-transformed, so exponentiate them
  mutate(trade_openness = exp(imports) + exp(exports)) %>% 
  mutate(gdp = exp(gdp)) %>% 
  select(iso3c, gdp, trade_openness)


## Merge back together
trade_gdp <- full_join(hold_dyadic %>% 
                          select(-imports_d_tradeflow, -imports_d_eite), 
                        imports_iso3_o, by = "iso3_o") %>% 
  # Get each country's total trade
  mutate(total_trade_o_tradeflow = exports_o_tradeflow + imports_o_tradeflow,
         total_trade_o_eite = exports_o_eite + imports_o_eite) %>% 
  # Now normalize trade based on GDP levels
  full_join(., gdp_levels, by = c("iso3_o" = "iso3c")) %>% 
  mutate(dyadic_tradeflow_normalized_gdp = dyadic_tradeflow/gdp,
         dyadic_eite_normalized_gdp = dyadic_trade_eite/gdp) %>% 
  select(iso3_o, iso3_d, 
         starts_with("dyadic"),
         starts_with("total")) 

## Make spatial weights with trade normalized to GDP 
spatial_weights_gdp <- left_join(trade_gdp, partners_targets_cs, by = "iso3_d") %>% 
  # Spatial weight in two steps
  rowwise() %>% 
  mutate(
    # -- Percentage targets
    multiply_percentage_tradeflow = ndc1_percentage_average_d * dyadic_tradeflow_normalized_gdp,
    multiply_percentage_eite = ndc1_percentage_average_d * dyadic_eite_normalized_gdp) %>% 
  group_by(iso3_o) %>% 
  summarize(
    # -- Percentage targets
    tradeflow_percentage = sum(multiply_percentage_tradeflow, na.rm = T),
    eite_percentage = sum(multiply_percentage_eite, na.rm = T)) %>% 
  ungroup()



#### Models for non-row-standardized spatial weights in *trade* ####

#' Order: 
#' -- Now that you've made the spatial weight for trade normalized by GDP:
#' -- Check correlation between both variables
#' -- Re-run baseline model for reference
#' -- Run model with new spatial weight (g1_a, g1_b)
#' -- Run model with EITE spatial weight (e1_a, e1_b)
#' -- Then run the interaction model (i1_a, i1_b)
#' -- Then for the interaction, show common support
#' -- Then for the interaction, show the fitted values
#' -- And make a table 

## Correlations between row-standardized and GDP-normalized spatial weight
data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, 
                     tradeflow_percentage_gdp = tradeflow_percentage, 
                     eite_percentage_gdp = eite_percentage),
            by = c("iso3c" = "iso3_o")) %>% 
  select(paris_target_safe, glasgow_target_safe, 
         tradeflow_percentage, tradeflow_percentage_gdp,
         eite_percentage, eite_percentage_gdp) %>% 
  cor(., use = "pairwise.complete") %>% 
  as.data.frame() %>% 
  select(ends_with("gdp")) 
#' correlation between total trade measures: 0.626
#' correlation between eite trade measures: 0.658

data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, 
                     tradeflow_percentage_gdp = tradeflow_percentage, 
                     eite_percentage_gdp = eite_percentage),
            by = c("iso3c" = "iso3_o")) %>% 
  select(paris_target_safe, glasgow_target_safe, 
         tradeflow_percentage, tradeflow_percentage_gdp,
         eite_percentage, eite_percentage_gdp) %>% 
  cor(., use = "complete") %>% 
  as.data.frame() %>% 
  select(ends_with("gdp"))
#' correlation between total trade measures: 0.836
#' correlation between eite trade measures: 0.772

## Run models

## "t1_a" and "t1_b" are the baseline models from the main results

## No controls 
t1_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage,
     data = .)

## With controls
t1_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)


## Now swap the GDP-normalized spatial weight
g1_a <- data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, tradeflow_percentage_gdp = tradeflow_percentage),
            by = c("iso3c" = "iso3_o")) %>% 
  mutate(tradeflow_percentage_gdp = tradeflow_percentage_gdp*1000) %>%
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage_gdp,
     data = .)

g1_b <- data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, tradeflow_percentage_gdp = tradeflow_percentage),
            by = c("iso3c" = "iso3_o")) %>%
  mutate(tradeflow_percentage_gdp = tradeflow_percentage_gdp*1000) %>%
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage_gdp +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)


## EITE trade
e1_a <- data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, eite_percentage_gdp = eite_percentage),
            by = c("iso3c" = "iso3_o")) %>% 
  mutate(eite_percentage_gdp = eite_percentage_gdp) %>%
  lm(glasgow_target_safe ~ paris_target_safe + eite_percentage_gdp,
     data = .)

e1_b <- data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, eite_percentage_gdp = eite_percentage),
            by = c("iso3c" = "iso3_o")) %>% 
  mutate(eite_percentage_gdp = eite_percentage_gdp) %>%
  lm(glasgow_target_safe ~ paris_target_safe + eite_percentage_gdp +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)
  
  
## Now interact trade openness with the row-standardized spatial weight
i1_a <- data_for_models %>%
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage*trade_openness,
     data = .)

i1_b <- data_for_models %>% 
  full_join(., spatial_weights_gdp %>% 
              select(iso3_o, tradeflow_percentage_gdp = tradeflow_percentage),
            by = c("iso3c" = "iso3_o")) %>%
  mutate(tradeflow_percentage_gdp = tradeflow_percentage_gdp*1000) %>%
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage*trade_openness +
       industry + renewables_share + fossil_rents,
     data = .)


## Common support between the row-standardized spatial weight and trade openness?

## Evaluate the model at the 17th and 83rd percentiles
quantile(data_for_models$trade_openness, probs = c(.25, .75)) # 3.99, 4.63
exp(quantile(data_for_models$trade_openness, probs = c(.25, .75))) # 54.3%, 103%

## For Low-Med-High; make boxplot to look for common support
quantile(data_for_models$trade_openness, probs = c(0.33, 0.67)) # 4.09, 4.56


## Boxplot for common support
common_support_trade <- data_for_models %>% 
  filter(!is.na(paris_target_safe), !is.na(glasgow_target_safe)) %>% 
  # Cut into terciles at 33rd and 67th (4.09 and 4.56)
  mutate(trade_openness_lmh = case_when(trade_openness < 4.09 ~ "Low",
                                        trade_openness > 4.09 & trade_openness < 4.56 ~ "Medium",
                                        trade_openness > 4.56 ~ "High")) %>% 
  mutate(trade_openness_lmh = fct_relevel(trade_openness_lmh, c("Low", "Medium", "High"))) %>% 
  ggplot(., aes(x = tradeflow_percentage, y = trade_openness_lmh, fill = trade_openness_lmh)) +
  labs(x = "Row-standardized spatial weight for trade", 
       y = "Trade openness (terciles)",
       fill = "Trade openness terciles") +
  scale_fill_manual(values = c("#F8766D", "#00BA38", "#619CFF")) +
  theme_classic() +
  theme(legend.position = "bottom") +
  geom_boxplot() +
  geom_point(shape = 1)

## Save 
ggsave(plot = common_support_trade, 
       filename = "plots/common-support-trade.pdf", 
       units = "in", height = 6, width = 8)


## Plot fitted values
#' Flat as high levels of openness, positive at low levels of openness
yhat_trade <- predictions(i1_a, 
            newdata = datagrid(tradeflow_percentage = seq(-47, 11, 3),
                               trade_openness = c(3.99, 4.63))) %>% 
  as.data.frame() %>% 
  select(estimate, starts_with("conf"), tradeflow_percentage, trade_openness) %>% 
  mutate(trade_openness = factor(case_when(trade_openness == 3.99 ~ "54% (25th percentile)", # Around the 17th percentile
                                           trade_openness == 4.63 ~ "103% (75th percentile)")), # Around the 90th percentile
         trade_openness = fct_relevel(trade_openness, "54% (25th percentile)", after = 0L)) %>%
  ggplot(., aes(tradeflow_percentage, estimate, fill = trade_openness)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(fill = "Trade openness", linetype = "Trade openness", colour = "Trade openness", 
       x = "Row-standardized spatial weight for trade", 
       y = "Glasgow target") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.5) +
  geom_line(aes(color = trade_openness, linetype = trade_openness), 
            alpha = 1, linewidth = 1.25)

## Save
ggsave(plot = yhat_trade, filename = "plots/yhat-trade.pdf", units = "in", height = 6, width = 8)


## Combine into a table

## Coefficient names
coefficients_renamed = c(
  "tradeflow_percentage" = "Trade Paris (row-standardized)",
  "tradeflow_percentage_gdp" = "Trade Paris (normalized by GDP)",
  "eite_percentage" = "Trade (EITE) Paris",
  "green_percentage" = "Trade (Green) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "eite_percentage_gdp" = "Trade (EITE) Paris (normalized by GDP)",
  "io_percentage" = "IO Paris (row-standardized)",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "paris_target_safe" = "Paris target",
  "io_percentage_non_rowstd" = "IO Paris (non--row-standardized)",
  "state_io_memberships" = "IO memberships",
  "(Intercept)" = "(Intercept)")

## Statistics to report
gm <- tibble::tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Observations", 0,
  "r.squared", "R2", 2)

## Combine in table
table_normalized_trade <- modelsummary(models = list(t1_a, t1_b, g1_a, g1_b, e1_a, e1_b, i1_a, i1_b),
             output = "latex_tabular",
             fmt = 2,
             coef_map = coefficients_renamed,
             coef_omit = c("industry|renewables_share|fossil_rents"),
             gof_map = gm,
             notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
#' -> paste into manuscript
#' Re-order rows (variables)
#' Add line space to controls
#' Bold the spatial matrix term
#' Add row with controls N/Y
#' Math the R2 term
#' Check caption


#### Models for non-row-standardized spatial weights in *IOs* ####

#' ToC
#' -- Make non-row-standardized spatial weight for IOs
#' -- Check correlations
#' -- Then run model
#' -- Then run the interaction
#' -- Then look for common support and show the fitted values
#' -- Then interact trade and IOs spatial weights; their common support; their fitted values

## Load IO membership data
io_membership_cs <- read_csv("input/io_membership_cs.csv") 

## Make non--row-standardized spatial weight
io_non_row_standardized <- full_imputed_averaged %>%
  # Subset
  select(iso3_o, iso3_d, sum_common_io) %>%
  # Group by dyad
  group_by(iso3_o, iso3_d) %>%
  # Average over years
  summarize(sum_common_io = mean(sum_common_io)) %>% 
  # Merge in targets
  left_join(., partners_targets_cs, by = "iso3_d") %>% 
  # Spatial weight in two steps
  rowwise() %>% 
  mutate(multiply_percentage_io = ndc1_percentage_average_d * sum_common_io) %>% 
  group_by(iso3_o) %>% 
  summarize(io_percentage_nonrowstd = sum(multiply_percentage_io, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(io_percentage_nonrowstd = io_percentage_nonrowstd/1000) %>% 
  rename(iso3c = iso3_o)

## Check correlation
full_join(io_non_row_standardized, data_for_models, by = "iso3c") %>% 
  select(paris_target_safe, glasgow_target_safe, io_percentage_nonrowstd, io_percentage) %>% 
  cor(., use = "pairwise.complete") %>% 
  as.data.frame() %>% 
  select(io_percentage) # 0.196

full_join(io_non_row_standardized, data_for_models, by = "iso3c") %>% 
  select(paris_target_safe, glasgow_target_safe, io_percentage_nonrowstd, io_percentage) %>% 
  cor(., use = "complete") %>% 
  as.data.frame() %>% 
  select(io_percentage) # 0.219

  


## Baseline models from table 2
t2_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage,
     data = .)
t2_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)

## Then remove the row-standardization by mutliplying the row-standardized spatial weight by state count of IO memberships
m2_a <- data_for_models %>%
  full_join(., io_non_row_standardized, by = "iso3c") %>%
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage_nonrowstd,
     data = .)

m2_b <- data_for_models %>% 
  full_join(., io_non_row_standardized, by = "iso3c") %>%
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage_nonrowstd +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)

## Interaction model
i2_a <- data_for_models %>% 
  full_join(io_membership_cs, by = "iso3c") %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage*state_io_memberships,
     data = .)

i2_b <- data_for_models %>% 
  full_join(io_membership_cs, by = "iso3c") %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage*state_io_memberships +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)


## Common support between the row-standardized spatial weight for IOs and IO memberships?

## Evaluate model at the 17th and 83rd percentiles 
quantile(data_for_models$io_percentage, probs = c(0.25, 0.75)) # -34.2, -29.3
data_for_models %>% 
  full_join(io_membership_cs, by = "iso3c") %>% 
  pull(state_io_memberships) %>% 
  quantile(., probs = c(0, .17, .25, .5, .75, .83, 1), na.rm = T) 
# 0 (min), 57.6 (17th), 62.5 (25th), 81 (50th), 103.5 (75th), 115.7 (83rd), 189 (max)


## Make boxplot with terciles to look at common support
common_support_ios <- data_for_models %>% 
  full_join(io_membership_cs, by = "iso3c") %>% 
  filter(!is.na(paris_target_safe), !is.na(glasgow_target_safe)) %>% 
  # Cut into terciles at 33rd and 67th (70 and 96)
  mutate(io_memberships_lmh = case_when(state_io_memberships < 70.1 ~ "Low",
                                       state_io_memberships > 70 & state_io_memberships < 96.1 ~ "Medium",
                                       state_io_memberships > 96 ~ "High")) %>% 
  mutate(io_memberships_lmh = fct_relevel(io_memberships_lmh, c("Low", "Medium", "High"))) %>% 
  ggplot(., aes(x = io_percentage, y = io_memberships_lmh, fill = io_memberships_lmh)) +
  labs(x = "Row-standardized spatial weight for IOs", 
       y = "Count of IO memberships (terciles)",
       fill = "Count of IO memberships terciles") +
  scale_fill_manual(values = c("#F8766D", "#00BA38", "#619CFF")) +
  theme_classic() +
  theme(legend.position = "bottom") +
  geom_boxplot() +
  geom_point(shape = 1)

## Save
ggsave(plot = common_support_ios, 
       filename = "plots/common-support-ios.pdf", units = "in", height = 6, width = 8)


## Plot fitted values 
#' Flat as high levels of membership, positive at low levels of membership
yhat_ios <- predictions(i2_a,
            newdata = datagrid(io_percentage = seq(-35, -23, 1),
                               state_io_memberships = c(63, 104))) %>%
  as.data.frame() %>%
  select(estimate, starts_with("conf"), io_percentage, state_io_memberships) %>%
  mutate(state_io_memberships = factor(case_when(state_io_memberships == 63 ~ "63 (25th percentile)", # Around the 17th percentile
                                                 state_io_memberships == 104 ~ "104 (75th percentile)")), # Around the 83th percentile
         state_io_memberships = fct_relevel(state_io_memberships, "63 (25th percentile)", after = 0L)) %>%  
  ggplot(., aes(io_percentage, estimate, fill = state_io_memberships)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(fill = "Count of IO memberships", 
       linetype = "Count of IO memberships",
       colour = "Count of IO memberships",
       x = "Row-standardized spatial weight for IOs",
       y = "Glasgow target") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.5) +
  geom_line(aes(color = state_io_memberships, linetype = state_io_memberships),
            alpha = 1, linewidth = 1.25)
## Save
ggsave(plot = yhat_ios, filename = "plots/yhat-ios.pdf", units = "in", height = 6, width = 8)


## Rename some coefficients
coefficients_renamed = c(
  "tradeflow_percentage" = "Trade Paris (row-standardized)",
  "tradeflow_percentage_gdp" = "Trade Paris (normalized by GDP)",
  "eite_percentage" = "Trade (EITE) Paris",
  "green_percentage" = "Trade (Green) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "io_percentage" = "IOs Paris (row-standardized)",
  "io_percentage_nonrowstd" = "IOs Paris (non--row-standardized)",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "state_io_memberships" = "IO memberships",
  "paris_target_safe" = "Paris target",
  "(Intercept)" = "(Intercept)")

## Combine in table
table_nonrowstd_io <- modelsummary(models = list(t2_a, t2_b, m2_a, m2_b, i2_a, i2_b),
             output = "latex_tabular",
             fmt = 2,
             coef_rename = coefficients_renamed,
             coef_omit = c("trade|industry|renewables_share|fossil_rents"),
             gof_map = gm,
             notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table_nonrowstd_io
#' -> paste into manuscript
#' Re-order rows (variables)
#' Add line space to controls
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption

## Where does the effect flatten out
predictions(i2_a,
            newdata = datagrid(io_percentage = seq(-35, -23, 1),
                               state_io_memberships = seq(125, 150, by = 1))) %>%
  as.data.frame() %>%
  select(estimate, io_percentage, state_io_memberships) %>% 
  filter(io_percentage %in% c(-35, -23)) %>% 
  arrange(state_io_memberships) %>% 
  pivot_wider(values_from = estimate, names_from = io_percentage, names_prefix = "io_perc") %>% 
  # Get difference in fitted values between low and high values of X, conditional on D
  mutate(diff = `io_perc-35` - `io_perc-23`) %>% 
  # Keep the point closest to zero
  filter(diff == min(abs(diff))) # 146 memberships

## What percentile is this?
ecdf(io_membership_cs$state_io_memberships)(146) # 93rd percentile...



#### Interact trade and IOs ####

## Additive models from the paper
t2_c <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage,
     data = .)
t2_d <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + tradeflow_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)

## Interaction models
trade_by_ios_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage*io_percentage,
     data = .)

trade_by_ios_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage*io_percentage +
       trade_openness + industry + renewables_share + fossil_rents,
     data = .)



## For Low-Med-High; make boxplot to look for common support

## Make boxplot to look for common support
quantile(data_for_models$tradeflow_percentage, probs = seq(0, 1, length = 7)) 
quantile(data_for_models$io_percentage, probs = seq(0, 1, length = 7)) 

## Make box plot
common_support_trade_ios <- data_for_models %>% 
  filter(!is.na(paris_target_safe), !is.na(glasgow_target_safe)) %>% 
  # Cut into terciles at 33rd and 67th
  mutate(tradeflow_percentage_lmh = case_when(tradeflow_percentage < -21.9 ~ "Low",
                                              tradeflow_percentage > -21.9 & tradeflow_percentage < -8.04 ~ "Medium",
                                              tradeflow_percentage > -8.04 ~ "High")) %>% 
  mutate(tradeflow_percentage_lmh = fct_relevel(tradeflow_percentage_lmh, c("Low", "Medium", "High"))) %>% 
  ggplot(., aes(x = io_percentage, y = tradeflow_percentage_lmh, fill = tradeflow_percentage_lmh)) +
  labs(x = "Spatial weight for IOs",
       y = "Spatial weight for trade (terciles)",
       # caption = "All spatial weights row-standardized",
       fill = "Spatial weight for trade terciles") +
  theme_classic() +
  theme(legend.position = "bottom") +
  geom_boxplot() +
  geom_point(shape = 1)

## Save
ggsave(plot = common_support_trade_ios, 
       filename = "plots/common-support-trade-ios.pdf", units = "in", height = 6, width = 8)


## Get 25th and 75th percentiles
quantile(data_for_models$io_percentage, probs = c(0.25, .75)) # -34, -29
quantile(data_for_models$tradeflow_percentage, probs = c(0.25, .75)) # -26, -5


## Fitted values, with trade on x-axis
yhat_trade_by_ios <- predictions(trade_by_ios_a, 
            newdata = datagrid(io_percentage = seq(-35, -24, by = 1),
                               tradeflow_percentage = seq(-47, 11, 1))) %>% 
  as.data.frame() %>% 
  select(estimate, starts_with("conf"), io_percentage, tradeflow_percentage) %>% 
  filter(io_percentage %in% c(-34, -29)) %>% 
  mutate(io_percentage = case_when(io_percentage == -34 ~ "-34% (25th percentile)",
                                   io_percentage == -29 ~ "-29% (75th percentile)"),
         io_percentage = fct_relevel(io_percentage, "-34% (25th percentile)", after = 0L)) %>% 
  ggplot(., aes(tradeflow_percentage, estimate, fill = io_percentage)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(fill = "Spatial weight for IOs", 
       linetype = "Spatial weight for IOs", 
       colour = "Spatial weight for IOs", 
       x = "Row-standardized spatial weight for trade", 
       y = "Glasgow target") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.5) +
  geom_line(aes(color = io_percentage, linetype = io_percentage), 
            alpha = 1, linewidth = 1.25)
## Save
ggsave(plot = yhat_trade_by_ios, filename = "plots/yhat-trade-by-ios.pdf", units = "in", height = 6, width = 8)


# Fitted values, with IOs on x-axis
predictions(trade_by_ios_a, 
            newdata = datagrid(io_percentage = seq(-35, -24, by = 1),
                               tradeflow_percentage = seq(-47, 11, 1))) %>% 
  as.data.frame() %>% 
  select(estimate, starts_with("conf"), io_percentage, tradeflow_percentage) %>% 
  filter(tradeflow_percentage %in% c(-26, -5)) %>% 
  mutate(tradeflow_percentage = case_when(tradeflow_percentage == -26 ~ "-26% (25th percentile)",
                                          tradeflow_percentage == -5 ~ "-5% (75th percentile)"),
         tradeflow_percentage = fct_relevel(tradeflow_percentage, "-26% (25th percentile)", after = 0L)) %>% 
  ggplot(., aes(io_percentage, estimate, fill = tradeflow_percentage)) +
  theme_classic() +
  theme(legend.position = "bottom") +
  labs(fill = "Spatial weight for trade", 
       linetype = "Spatial weight for trade", 
       colour = "Spatial weight for trade", 
       x = "Row-standardized spatial weight for IOs", 
       y = "Glasgow target") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high),
              alpha = 0.5) +
  geom_line(aes(color = tradeflow_percentage, linetype = tradeflow_percentage), 
            alpha = 1, linewidth = 1.25)
## Save
# ggsave(plot = last_plot(), filename = "plots/yhat-ios-by-trade.pdf", units = "in", height = 6, width = 8)




## Combine in table
table_trade_by_ios <- modelsummary(models = list(t2_c, t2_d, 
                           trade_by_ios_a, trade_by_ios_b),
             output = "latex_tabular",
             fmt = 2,
             coef_rename = coefficients_renamed,
             coef_omit = c("industry|renewables_share|fossil_rents"),
             gof_map = gm,
             notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table_trade_by_ios
#' -> paste into manuscript
#' Re-order rows (variables)
#' Add line space to controls
#' Bold the spatial matrix term
#' Math the R2 term
#' Check caption





#### Regions ####

## Control for unobserved regional heterogeneity
r1_a <- data_for_models %>%  
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage + 
       as.factor(region), # East Asia is the reference category 
     data = .)

r1_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage + 
       trade_openness + industry + renewables_share + fossil_rents +
       as.factor(region), # East Asia is the reference category 
     data = .)

r2_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + 
       as.factor(region), # East Asia is the reference category 
     data = .)

r2_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + io_percentage + 
       trade_openness + industry + renewables_share + fossil_rents +
       as.factor(region), # East Asia is the reference category 
     data = .)

r3_a <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage + io_percentage + 
       as.factor(region), # East Asia is the reference category 
     data = .)

r3_b <- data_for_models %>% 
  lm(glasgow_target_safe ~ paris_target_safe + tradeflow_percentage + io_percentage + 
       trade_openness + industry + renewables_share + fossil_rents +
       as.factor(region), # East Asia is the reference category 
     data = .)

## Statistics to report
gm <- tibble::tribble(
  ~raw, ~clean, ~fmt,
  "nobs", "Observations", 0,
  "r.squared", "R2", 2)

coefficients_renamed = c(
  "tradeflow_percentage" = "Trade Paris (row-standardized)",
  "tradeflow_percentage_gdp" = "Trade Paris (normalized by GDP)",
  "eite_percentage" = "Trade (EITE) Paris",
  "green_percentage" = "Trade (Green) Paris",
  "competition_tradeflow_percentage" = "Trade (competition) Paris",
  "competition_eite_percentage" = "Trade (competition, EITE) Paris",
  "io_percentage" = "IO Paris (row-standardized)",
  "trade_openness" = "Trade openness",
  "renewables_share" = "Renewable electricity",
  "industry" = "Industry",
  "fossil_rents" = "Fossil rents",
  "as.factor(region)Europe & Central Asia" = "Europe and Central Asia",
  "as.factor(region)Latin America & Caribbean" = "Latin America and Caribbean",
  "as.factor(region)Middle East & North Africa" = "Middle East and North Africa",
  "as.factor(region)North America" = "North America",
  "as.factor(region)South Asia" = "South Asia",
  "as.factor(region)Sub-Saharan Africa" = "Sub-Saharan Africa",
  "paris_target_safe" = "Paris target",
  "io_percentage_non_rowstd" = "IO Paris (non--row-standardized)",
  "state_io_memberships" = "IO memberships",
  "(Intercept)" = "(Intercept)")

table_regions <- modelsummary(models = list(r1_a, r1_b, r2_a, r2_b, r3_a, r3_b),
             output = "latex_tabular",
             fmt = 3,
             coef_rename = coefficients_renamed,
             coef_omit = c("trade_openness|industry|renewables_share|fossil_rents"),
             gof_map = gm,
             statistic = 'p.value',
             notes = list("Outcome is Glasgow mitigation target"),
             stars = c("+" = 0.1, "*" = 0.05, "**" = 0.01))
table_regions


